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Abstract 

The  Generalized  Empirical  Interpolation  Method  (GEIM)  is  an  extension  first  presented  by  Maday  and  Mula  in  [1]  in 
2013  of  the  classical  empirical  interpolation  method  (presented  in  2004  by  Barrault,  Maday,  Nguyen  and  Patera  in  [2]) 
where  the  evaluation  at  interpolating  points  is  replaced  by  the  more  practical  evaluation  at  interpolating  continuous 
linear  functionals  on  a  class  of  Banach  spaces.  As  outlined  in  [1],  this  allows  to  relax  the  continuity  constraint  in  the 
target  functions  and  expand  both  the  application  domain  and  the  stability  of  the  approach.  In  this  paper,  we  present 
a  thorough  analysis  of  the  concept  of  stability  condition  of  the  generalized  interpolant  (the  Lebesgue  constant)  by 
relating  it  to  an  inf-sup  problem  in  the  case  of  Hilbert  spaces.  In  the  second  part  of  the  paper,  it  will  be  explained  how 
GEIM  can  be  employed  to  monitor  in  real  time  physical  experiments  by  providing  an  online  accurate  approximation  of 
the  phenomenon  that  is  computed  by  combining  the  acquisition  of  a  minimal  number,  optimally  placed,  measurements 
from  the  processes  with  their  mathematical  models  (parameter-dependent  PDEs).  This  idea  is  illustrated  through  a 
parameter  dependent  Stokes  problem  in  which  it  is  shown  that  the  pressure  and  velocity  fields  can  efficiently  be 
reconstructed  with  a  relatively  low-dimensional  interpolation  space. 

Keywords:  empirical  interpolation;  generalized  empirical  interpolation;  reduced  basis;  model  order  reduction; 
stability;  Stokes  equations 


Introduction 

Let  A  be  a  Banach  space  of  functions  defined  over  a  domain  Q  c  (or  C”^).  Let  (A„)„(=n.,  A„  c  A,  be  a  family 
of  finite  dimensional  spaces,  dim  A„  =  n,  and  let  be  an  associated  family  of  sets  of  points:  S„  -  {v")'Lj,  with 

x"  e  Q.  The  problem  of  interpolating  any  function  /  e  A  has  traditionally  been  stated  as: 

”Eind  /„  e  A„  such  that  f„(x")  =  /(x”),  Vi  e  {1, ... , «)”,  (1) 

where  we  note  that  it  is  implicitly  required  that  A  is  a  Banach  space  of  continuous  functions.  The  most  usual  ap¬ 
proximation  in  this  sense  is  the  Lagrangian  interpolation,  where  the  interpolating  spaces  X„  are  of  polynomial  nature 
(spanned  by  plain  polynomials,  rational  functions,  Eourier  series...)  and  the  question  on  how  to  appropriately  select 
the  interpolating  points  in  this  case  has  broadly  been  explored.  Although  there  exists  still  nowadays  open  issues  on 
Lagrangian  interpolation  (see,  e.g.  [3]),  it  is  also  interesting  to  look  for  extensions  of  this  procedure  in  which  the  in¬ 
terpolating  spaces  X„  are  not  necessarily  of  polynomial  nature.  The  search  for  new  interpolating  spaces  A„  is  therefore 
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linked  with  the  question  on  how  to  optimally  select  the  interpolating  points  in  this  case  and  how  to  obtain  a  process 
that  is  at  least  stable  and  close  to  the  best  approximation  in  some  sense. 

Although  several  procedures  have  been  explored  in  this  direction  (we  refer  to  [4],  [5]  and  also  to  the  kriging  studies 
in  the  stochastic  community  such  as  [6]),  of  particular  interest  for  the  present  work  is  the  Empirical  Interpolation 
Method  (EIM,  [2],  [7],  [8])  that  has  been  developed  in  the  broad  framework  where  the  functions  /  to  approximate 
belong  to  a  compact  set  F  of  continuous  functions  {X  -  C°(Q)).  The  structure  of  F  is  supposed  to  make  any  f  &  F 
be  approximable  by  finite  expansions  of  small  size.  This  is  quantified  by  the  Kolmogorov  «-width  dn{F,X)  of  F  in 
X  (see  definition  (2)  below)  whose  smallness  measures  the  extent  to  which  F  can  be  approximated  by  some  finite 
dimensional  space  X„  c  X  of  dimension  n.  Unfortunately,  in  general,  the  best  approximation  «-dimensional  space 
is  not  known  and,  in  this  context,  the  Empirical  Interpolation  Method  aims  to  build  a  family  of  interpolating  spaces 
X„  with  satisfactory  approximation  properties  together  with  sets  of  interpolating  points  S„  such  that  the  interpolation 
is  well  posed.  This  is  done  by  a  greedy  algorithm  on  both  the  interpolating  points  and  the  interpolating  selected 
functions  ipi  (see  [2]).  This  procedure  has  the  main  advantage  of  being  constructive,  i.e.  the  sequence  of  interpolating 
spaces  iXn)  and  interpolating  points  (5„)  are  hierarchically  defined  and  the  procedure  can  easily  be  implemented  by 
recursion. 

A  recent  extension  of  this  interpolation  process  consists  in  generalizing  the  evaluation  at  interpolating  points  by 
application  of  a  class  of  interpolating  continuous  linear  functionals  chosen  in  a  given  dictionary  E  c  £,{X).  This  gives 
rise  to  the  so-called  Generalized  Empirical  Interpolation  Method  (GEIM).  In  this  new  framework,  the  particular  case 
where  the  space  X  -  L}{Q)  was  first  studied  in  [1].  We  also  mention  the  preliminary  works  of  [9]  in  which  the  authors 
introduced  the  use  of  linear  functionals  in  EIM  in  a  finite  dimensional  framework.  In  the  present  paper,  we  will  start 
by  revisiting  the  foundations  of  the  theory  in  order  to  show  that  GEIM  holds  for  Banach  spaces  X  (Section  1).  The 
concept  of  stability  condition  (Lebesgue  constant,  A„)  of  the  generalized  interpolant  will  also  be  introduced. 

In  the  particular  case  where  A  is  a  Hilbert  space,  we  will  provide  an  interpretation  of  the  generalized  interpolant  of 
a  function  as  an  oblique  projection.  This  will  shed  some  light  in  the  understanding  of  GEIM  from  an  approximation 
theory  perspective  (Section  2.1).  This  point  of  view  will  be  the  key  to  show  that  the  Lebesgue  constant  is  related  to 
an  inf-sup  problem  (Section  2.2)  that  can  be  easily  computed  (Section  3).  The  derived  formula  can  be  seen  as  an 
extension  of  the  classical  formula  for  Lagrangian  interpolation  to  Hilbert  spaces.  It  will  also  be  shown  that  the  Greedy 
algorithm  aims  to  minimize  the  Lebesgue  constant  in  a  sense  that  will  be  made  precise  in  Section  2.3.  Eurthermore, 
the  inf-sup  formula  that  will  be  introduced  will  explicitly  show  that  there  exists  an  interaction  between  the  dictionary 
E  of  linear  functionals  and  the  Lebesgue  constant.  Although  it  has  so  far  not  been  possible  to  derive  a  general  theory 
about  the  impact  of  E  on  the  behavior  of  the  Lebesgue  constant,  we  present  in  Section  4  a  first  simple  example  in 
which  this  influence  is  analyzed  through  numerical  simulation. 

The  last  part  of  the  paper  (Section  5)  will  allow  to  present  some  more  elaborate  potential  applications  of  the 
method  with  respect  to  what  is  presented  in  [1].  In  particular,  we  will  explain  how  GEIM  can  be  used  to  build  a  tool 
for  the  real-time  monitoring  of  a  physical  or  industrial  process.  This  will  be  achieved  thanks  to  the  online  computation 
of  a  generalized  interpolant  that  will  approximate  the  phenomenon  under  consideration.  Its  derivation  will  combine 
measurements  collected  on  the  fly  from  the  process  itself  with  a  mathematical  model  (a  parameter  dependent  PDE) 
that  represents  the  physical  understanding  of  the  process.  It  will  also  be  explained  how  the  proposed  methodology 
can  be  helpful  for  the  minimization  of  the  number  of  sensors  required  to  reconstruct  the  held  variable  and  also  their 
optimal  selection  and  placement,  which  are  very  important  issues  in  engineering.  These  ideas  will  be  illustrated 
through  a  parameter  dependent  Stokes  problem  for  X  -  (//'(Q))  x  Lq(Q),  where  Lq(Q)  is  the  space  of  the  L^(Q) 
functions  with  zero  mean  over  Q. 

Taking  advantage  of  this  idea,  we  will  outline  in  the  conclusion  how  the  method  could  be  used  to  build  an  adaptive 
tool  for  the  supervision  of  experiments  that  could  distinguish  between  normal  and  accidental  conditions.  We  believe 
that  this  tool  could  help  in  taking  real-time  decisions  regarding  the  security  of  processes. 

1.  The  Generalized  Empirical  Interpolation  Method 

Let  A  be  a  Banach  space  of  functions  defined  over  a  domain  Q.  c  where  d  -  1,2, 3.  Its  norm  is  denoted  by 
IMU.  Let  F  he  a  compact  set  of  A.  With  M  being  some  given  large  number,  we  assume  that  the  dimension  of  the 
vectorial  space  spanned  by  F  (denoted  asT  -  span{E))  is  of  dimension  larger  than  M.  Our  goal  is  to  build  a  family 
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of  n-dimensional  subspaces  of  X  that  approximate  well  enough  any  element  of  F.  The  rationale  of  this  approach  is 
linked  to  the  notion  of  n-width  following  Kolmogorov  [10]: 

Definition  1.1.  Let  F  be  a  subset  of  some  Banach  space  X  and  ¥„  be  a  generic  n-dimensional  subspace  of  X.  The 
deviation  between  F  and  Y„  is 

E(F-,  Y„)  sup  inf  ||;c  -  y|U  ■ 

xeF  y<^Y„ 

The  Kolmogorov  n-width  ofF  in  X  is  given  by 

d„(F,X)  inf{£’(T';  Y„)  :  ¥„  a  n-dimensional  subspace  of  X] 

=  inf  sup  inf  ||x-y|U  ■  (2) 

Y„<lX 
dim  Y„=n 

The  smallness  of  the  n-width  of  F  thus  measures  to  what  extent  the  set  F  can  be  approximated  by  an  n-dimensional 
subspace  of  X.  Several  reasons  can  account  for  a  rapid  decrease  of  d„{F,X):  if  T  is  a  set  of  functions  defined  over 
a  domain,  we  can  refer  to  regularity,  or  even  to  analyticity,  of  these  functions  with  respect  to  the  domain  variable  (as 
analyzed  in  the  example  in  [10]).  Another  possibility  —  that  will  actually  be  used  in  our  numerical  application —  is 
when  F  -  {u(p, .),  p  6  D},  where  D  is  a  compact  set  of  and  u{p, .)  is  the  solution  of  a  PDE  parametrized  by  p. 
The  approximation  of  any  element  u(p, .)  e  F  hy  finite  expansions  is  a  classical  problem  addressed  by,  among  others, 
reduced  basis  methods  and  the  regularity  of  u  in  p  can  also  be  a  reason  for  having  a  small  n-width  as  the  results  of 
[11]  and  [12]  show. 

Finally,  let  us  also  assume  that  we  have  at  our  disposal  a  dictionary  of  linear  functionals  F  c  £.{X)  with  the 
following  properties: 

PI:  VcreF,|Mlx(.Y)  =  l. 

P2:  Unisolvence  property:  If  (,0  e  spanlF]  is  such  that  cr{(f)  -  0,  Vcr  e  S,  then  (f  -0. 

Given  this  setting,  GEIM  aims  at  building  M-dimensional  interpolating  spaces  Xm  spanned  by  suitably  chosen  func¬ 
tions  {ipi,<p2,  ■  ■  ■ ,  ^m]  of  F  together  with  sets  of  M  selected  linear  functionals  (cri,  cr2,. . . ,  ctm]  coming  from  E  such 
that  any  ip  e  F  k  well  approximated  by  its  generalized  interpolant  J7M[y5]  £  Xm  defined  by  the  following  interpolation 
property: 

Tp  e  X,  Jm[p\  £  Xm  such  that  crfjMip])  -  a-fp),  V/  =  1, . . . ,  M.  (3) 

Remark  1.2.  Since  only  some  elements  of  the  dictionary  E  are  going  to  be  selected,  note  that  E  consists  of  ’’candidate  ” 
functionals  and  only  the  selected  elements  will  actually  be  implemented. 

The  definition  of  GEIM  in  the  sense  of  (3)  raises  several  questions: 

•  is  there  an  optimal  selection  for  the  linear  functionals  cr,  within  the  dictionary  E  ? 

•  is  there  a  constructive  optimal  selection  for  the  functions  pt  6  FI 

•  given  a  set  of  linearly  independent  functions  {pi]ie[\,M]  and  a  set  of  continuous  linear  functionals  (cr,],g[i  does 

the  interpolant  exist  in  the  sense  of  (3)? 

•  is  the  interpolant  unique? 

•  under  what  hypothesis  can  we  expect  the  GEIM  approximation  to  converge  rapidly  to  pi 

In  what  follows,  we  provide  answers  to  these  questions  either  with  rigorous  proofs  or  with  numerical  evidences. 

The  construction  of  the  generalized  interpolation  spaces  Xm  and  the  selection  of  the  suitable  associated  linear 
functionals  is  recursively  performed  by  following  a  greedy  procedure  very  similar  to  the  one  of  the  classical  EIM.  The 
first  selected  function  is,  e.g., 

py  =  argsup||(,c||^, 

tpeF 
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that  defines  Xi  -  span{(^i).  The  first  interpolating  linear  functional  is 

0-1  =  argsup  \a-{ipx)\. 

creL 

The  interpolation  operator  SF\  ;  is  defined  such  that  (3)  is  true  for  M  -  i.e.  cri  -  crx{ip),  for  any 

i/j  6  /Y.  To  facilitate  the  practical  computation  of  the  generalized  interpolant,  we  express  it  in  terms  of 

<7l  =  - , 

which  will  be  the  basis  function  that  will  be  employed  for  Xi .  In  this  basis,  the  interpolant  reads 

J'i[V5]  =  0-1(9?)^  1,  SipeX. 

We  then  proceed  by  induction.  For  Mtnax  <  A1  an  upper  bound  fixed  a  priori,  assume  that,  for  a  given  \  <  M  < 
we  have  selected  a  set  of  functions  {931 , 932, ... ,  93m)  and  the  associated  basis  functions  {q\,q2,  ■  ■  ■  ^qm]  that  span  Xm, 
as  well  as  the  interpolating  linear  functionals  {cri,cr2, . . .  ,o"m).  The  generalized  interpolant  is  assumed  to  be  well 
defined  by  (3),  i.e., 

M 

j=i 

where  the  coefficients  j  -  1, . . . ,  M  are  given  by  the  interpolation  problem 

Find  {a^(93))^i  such  that; 

,  M 

U=i 

where  are  the  coefficients  of  the  M  x  M  matrix  :=  (crTo,))  .  We  now  define 

V90  e  F,  SMip)  =  Wp  -  JmMIU- 

At  the  M  +  1-th  stage  of  the  greedy  algorithm,  we  choose  ipM+\  such  that 


Pm+i  =  argsup  Em (95) 

ipeF 


(4) 


and 


ctm+i  =  argsup  |cr(93M+i  -  J'm[9^m+i])I- 

(T€£ 


(5) 


The  next  basis  function  is  then 


_  ‘Pm+i  -  STmIpm+i] 

(JM+1  ~  -  • 

0"M+i(i,OM+1  -  ,!7m[9^m+i]) 

We  finally  set  Xm+\  =  span{(^9,  1<  j  <  M  +  \]  -  spanj^^,  1  <  j  <  M  +  \].  The  interpolation  operator  J'm+i  '■  X  i-» 
^M+i  is  given  by 


M+l 

so  as  to  satisfy  (3).  The  coefficients  j  -  .  ,M  +  1,  are  therefore  given  by  the  interpolation  problem 


Find  such  that; 

M+l  ^ 

Z  =  o-i(<p),  Vi  =  1, . . . ,  M  -H  1, 

J=i 
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where  =  (cr,(fl,)) 

By  following  exactly  the  same  guidelines  as  in  [1]  where  the  particular  case  X  =  L^(Q)  was  addressed,  it  can 
be  proven  that,  in  the  general  case  where  ^  is  a  Banach  space,  the  generalized  interpolation  is  well-posed:  for  any 
1  <  M  <  Ai,  the  set  of  functions  {qj,  j  e  [1,  M])  is  linearly  independent  and  therefore  the  space  Xm  is  of  dimension 
M.  Furthermore,  the  matrix  is  lower  triangular  with  unity  diagonal  (hence  invertible)  with  off-diagonal  entries  in 
[-1,1]. 

Note  that  GEIM  reduces  to  EIM  if  /Y  c  C°(f2)  and  2  is  composed  of  Dirac  masses.  Also,  if  the  cardinality  #F  of  F 
is  finite,  then  the  Greedy  algorithm  is  exact  in  the  sense  that  F  c  X#f.  This  type  of  property  does  not  hold  in  traditional 
Lagrangian  interpolation  due  to  the  fact  that  the  interpolating  polynomial  spaces  are  used  to  interpolate  continuous 
functions  that  are  not  necessarily  of  polynomial  nature.  Einally,  note  also  that  the  approach  can  be  shortcut  if  the  basis 
functions  are  available,  in  which  case  the  interpolating  linear  functionals/points  are  the  only  output  of  GEIM/EIM. 

It  is  also  important  to  point  out  that  the  current  extension  of  EIM  presents  two  major  advantages:  first,  it  allows  the 
interpolation  of  functions  of  weaker  regularity  than  C*’(Q).  The  second  interest  is  related  to  the  potential  applications 
of  GEIM:  the  use  of  linear  functionals  can  model  in  a  more  faithful  manner  real  sensors  involved  in  physical  exper¬ 
iments  (indeed,  these  are  in  practice  no  point  evaluations  as  it  is  usually  supposed  but  rather  local  averages  of  some 
quantity  of  interest).  The  potentialities  of  these  two  aspects  will  be  illustrated  in  the  numerical  application  presented 
in  Section  5. 

We  now  state  a  first  result  about  the  interpolation  error  of  GEIM. 

Theorem  1.3  (Interpolation  error  on  a  Banach  space).  Vi/:  e  X,  the  interpolation  error  satisfies 


\\<f>  -  SF m[‘P\\\x  ^  +  ^m)  inf  W^' 


'PmWx^ 


where 


Am  II.!7mIIx(A’)  -  sup 

ipeX 


IIJ'mMIU 

ll^lk 


is  the  Lebesgue  constant  in  the  X  norm. 


Proof.  The  desired  result  easily  follows  since  for  any  ip  e  X  and  any  6  Xm  we  have 


< 

< 


\\[(p  -  tfu]  -  FTMif  -  iAmIIU 
lllf  -  SfM\\t{X)\\v  -  'FmWx 
(1  +  \VTm\\j:(X))\\v  -  </'mIU, 


which  yields  the  desired  inequality. 


(6) 

(7) 


□ 


The  last  term  in  the  right  hand  side  of  equation  (6)  is  known  as  the  best  fit  of  p  by  elements  in  the  space  Xm. 
However,  Xm  does  not  in  general  coincide  with  the  optimal  M-dimensional  space  in  the  sense  that  Xm  +  X°^\  with 

X°M  ^  arginf  E(F,Ym). 

YmcX 

dim{YM)=M 

This  raises  the  question  of  the  quality  of  the  finite  dimensional  subspaces  Xm  provided  by  the  Greedy  selection 
procedure.  It  has  been  proven  first  in  [13]  in  the  case  of  X  -  LF{Q.)  and  then  in  [14]  in  a  general  Banach  space  that 
the  interpolating  spaces  Xm  coming  from  the  Greedy  selection  procedure  of  GEIM  are  quite  optimal  and  that  the  lack 
of  optimality  comes  from  the  Lebesgue  constant.  The  main  results  are  the  following  (see  [14]): 

Theorem  1.4  (See  Corollary  3.13  of  [14]). 


i)  If  dM(F,X)  <  CqM  M  —  1,2, . . .  and  that  (1  -H  Am)  <  C^Mffor  any  M  —  1,2,...,  then,  for  all  f  >  1/2,  the 
interpolation  error  satisfies  for  any  ip  e  F  the  inequality  \\ip  —  fTMW\\\x  ^  where 


Cl 


max 


1 

2  2 


max 


M=1 


max  M“  ^  ^ 

...,2L2(f+/S)J+l 
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ii)  If  (Am)  is  a  monotonically  increasing  sequence  and  if  dM(F,X)  <  Coe 
then,  for  any  ip  €  F,  the  interpolation  error  can  be  bounded  as 


for  any  M  >  I  and  with  Co  >  1, 


W'P  -  ffM[‘P]\\x  ^ 


|4Co(l+Ai), 

I  V2Co(l  +  Am)^ 


1. 

ifM  >  2. 


As  a  consequence  of  this  result,  the  interpolation  error  of  GEIM  will  converge  if  the  Lebesgue  constant  ’’does  not 
increase  too  fast”  in  the  sense  that  it  allows  that  the  previous  upper  bounds  tend  to  zero  as  the  dimension  M  increases. 
By  following  the  same  lines  as  in  [1],  it  can  be  proven  that  when  A  is  a  Banach  space,  the  Lebesgue  constant  has  the 
exponential  upper-bound 

Am  <  2"-'  max  ||^,|U,  (8) 

which  implies  that  the  decay  of  dM(F,  X)  should  be  exponential  in  order  to  converge.  However,  the  behavior  of  (Am)m 
observed  in  numerical  applications  (see  Section  5)  is  rather  linear  and  leads  us  to  expect  that  the  upper  bound  of  (8)  is 
far  from  being  optimal  in  a  class  of  set  F  of  small  Kolmogorov  n-width. 


2.  Further  results  in  the  case  of  a  Hilhert  space 

In  this  section  A  is  a  Hilbert  space  of  functions  where  the  norm  IMU  is  induced  by  the  inner  product  (■,  ■)a'-  We  will 
see  that  in  this  case  the  generalized  interpolant  can  be  seen  as  an  oblique  projection.  It  will  also  be  proven  that  we  can 
derive  a  sharp  interpolation  error  bound  in  this  case.  An  explicit  (and  easily  computable)  formula  for  the  Lebesgue 
constant  will  also  be  obtained  and  this  formula  will  be  used  to  show  that  the  Greedy  algorithm  aims  at  minimizing  the 
Lebesgue  constant. 


2.1.  Interpretation  of  GEIM  as  an  oblique  projection 

Lor  1  <  j  <  M,  if  cTj  is  the  /^-linear  functional  selected  by  the  greedy  algorithm,  let  w,  be  its  Riesz  representation 
in  A,  i.e.  wj  is  such  that 

V/eA,  (rj(f)^(wjJ)x.  (9) 

It  follows  from  the  well  posedness  of  the  generalized  interpolation  that  {cri, . . .  ,0"^)  are  linearly  independent  and 
therefore  (wi, . . . ,  wm)  are  also  linearly  independent.  We  will  denote  by  Wm  the  M-dimensional  space 

Wm  =  span{wi,...,WM). 

Lor  any  /  6  A,  let  nn/j^[/]  be  the  orthogonal  projection  of  /  on  Wm,  i.e. 

|nw„[/]  e  Wm 

[(/- nn?„[/],w)A-  =  0,  Vw  e  Wm- 

With  these  notations,  we  can  provide  the  following  interpretation  of  the  generalized  interpolant  of  a  function  (see 
Ligure  1  for  a  schematic  representation): 

Lemma  2.1.  V/  e  A,  ffM[f]  is  an  oblique  projection  onto  the  space  Xm  orthogonal  to  the  space  Wm,  i-e. 


[jM[f]^XM 

|(J'm[/]-/,w).y=0,  Vwe  Wm. 

Proof.  Lor  any  /  e  A,  the  interpolation  property  reads  crj(f)  -  crjifTMiff),  for  1  <  j  <  M.  It  is  then  clear  that 
(w j,  f)x  —  (wj,  ffMifVtx  and  the  result  easily  follows  from  the  fact  that  {wi ,...,  wm)  are  a  basis  of  Wm-  □ 

A  direct  consequence  of  Lemma  2. 1  is  the  following  result: 
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/ 


Figure  1:  Interpretation  of  J7m[/]  as  an  oblique  projection. 


Corollary  2.2,  In  the  particular  case  where 

2  =  {o-  6  XW,  s.t.  V/  e  span{F)-",  cr(/)  =  o) , 

then  Wm  —  ond  the  resulting  generalized  interpolant  is  the  orthogonal  projection  of  f  onto  the  space  Xm,  i-e. 

JmU]  -  An- 

Proof.  First  of  all,  note  that,  in  this  setting,  there  is  a  bijective  mapping  between  2  and  span{F)  (because  the  Riesz 
representation  of  any  cr  e  2  is  an  element  of  span{F)  and  vice-versa).  Now,  from  the  argmax  definition  of  ctj.  in  the 
greedy  algorithm,  the  Riesz  representation  of  crj.  is  the  function  Wk  -  <fk-  Ak-\{<fk\  for  k  >2  and  wi  -  ip\\fk  -  1. 
The  interpolation  property  CTkif  -  Ahiin)  —  0  implies  in  this  case  that  {wk,f  -  fTk-\{fY)x  for  any  k  e  {1, . . .  ,M). 
But  since  the  family  {wi, . . . ,  wm}  is  a  basis  of  Xm  in  this  particular  case,  it  follows  that  (/  -  fTAtilf],  w)x  -  0  for  all 
w  6  Xm.  C 

Remark  2.3.  The  case  E  =  {cr  e  X  (<T) ,  s.t.  V/  6  span{F')'‘-,  cr(f)  —  0)  is  a  theoretical  situation  that  does  not  usually 
hold  in  practical  applications.  Corollary  2.2  is  however  a  first  step  towards  the  theoretical  understanding  of  the 
impact  of  the  dictionary  E  on  the  interpolation  procedure. 

From  Lemma  2.1,  note  that  fTmlf]  can  also  be  seen  as  a  particular  Petrov-Galerkin  approximation  of  the  function 
/  in  the  case  where  the  approximation  space  is  Xm  and  the  trial  space  is  Wm-  Indeed,  the  search  for  the  generalized 
interpolant  can  be  stated  as 

jOiven  f  e  X,  find  JmU]  6  Xm  such  that 
\iJM[n,^)x  ^  UwU,  -iweWM. 

This  formulation  leads  to  the  classical  error  estimation; 


II/-J'm[/]IU< 


inf  II/- 


(13) 


where  /m  is  the  inf-sup  constant 


Pm  ;=  inf 

xeXm 


ix,  w)x 

sup - . 

weWmWxWxWwWx 


(14) 


It  will  be  proven  in  the  next  section  that  the  parameter  IIPm  is,  in  fact,  equal  to  the  Lebesgue  constant  A^.  We  will 
also  see  that  the  error  bound  provided  in  relation  (13)  is  slightly  suboptimal  due  to  the  presence  of  the  coefficient  1 
before  the  parameter  \  j Pm- 


2.2.  Interpolation  error 

The  interpretation  of  the  generalized  interpolant  as  an  oblique  projection  is  useful  to  derive  the  following  result 
about  the  interpolation  error: 
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Theorem  2.4  (Interpolation  error  on  a  Hilbert  space).  Sip  e  X,  the  interpolation  error  satisfies  the  sharp  upper  bound: 


\\V-JmW\\\x^^m  inf 


(15) 


where  Am  \\Jm\\z(X)  =  sup 


\\JmW]\\x 

IMU 


is  the  Lebesgue  constant  in  the  X  norm.  Furthermore,  Am  =  — ,  where 

Pm 


Pm  inf  sup - . 

{W'‘~,  fi^)x 

Proof.  LetVM:=  inf  sup  - — — — — It  immediately  follows  that 

\\x\W  \\x 


(16) 


,f)x 


Vw  eWM,  vm\\w  IU<  sup 

W  \\x 


Furthermore,  for  any  p  &  X,\i  follows  from  Lemma  2.1  that  ip  -  £  ^m-  Then: 

II  rr  r  nil  i^-JMW],>P^)x 

vmW^- JM{ip\\\x  ^  sup  - — — - . 

\W  \\X 

Besides,  for  any  e  Xm  and  any  6  X^: 

-  JmVp\,  >P^)x  ^^)x  ■ 

The  Cauchy-Schwarz  inequality  applied  to  (17)  combined  with  relation  (18)  yields: 


(17) 


(18) 


Vm\W-  jM[‘f]\\x  <  ttxfWp-fiWx- 

tpeXM 

Next,  it  can  be  proven  (see  the  proposition  of  Appendix  A.l)  that  vm  -  Pm,  which  yields  the  inequality 

\\(p  -  JmW]\\x  <  inf  \\(p  -  tfWx- 

pi^  (peXM 

The  end  of  the  proof  consists  in  showing  that 


(19) 


(20) 


1  *  IIJ'mMIU 

—  =  =  sup—— - 

Pm  ipeX  wWx 

This  is  done  by  noting  first  of  all  that  formula  (16)  implies  that 


€  A,  PM\\jMm\x  <  sup  <  ll^ll^, 

weWu  Iln’II-f 


(21) 


where  we  have  used  the  fact  that  {fJMVp\,w)x  -  i.'p,  w)x  for  all  w  e  Wm  and  the  Cauchy-Schwarz  inequality.  There¬ 
fore, 

'ip  e  A,  IIJ-mMIU  < 

Pm 


which  yields 


Am 


(22) 
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Let  us  now  denote  by  ijj  an  element  of  Xm  with  norm  ||i^|U  =  1  such  that 


sup 


(i//,  Wm)x 

IIw^mIU 


-  Pm- 


Using  the  notation  introduced  in  (10),  we  denote  by  T\-Wm  orthogonal  projection  oitj/  owWm-  Similarly,  [i/'] 

is  the  orthogonal  projection  of  i/r  on  W-^  so  that  -  0 

^  +  nvr^[i/']. 


Note  that 

SFm  [njv„[iA]]  =  <A- 

Indeed,  since  i/r  6  Xm,  we  have  that 

lA  =  PTm  [iA]  =  PTm  [nn/„[iA]]  +  J^M  [rityx  [^]j  =  J'm  [niy„[iA]] 
because  J'm  [niv^  [lAl]  =  0  in  vertue  of  the  interpolation  property  given  in  relation  (11).  In  addition  to  this, 

(lA,  wm)x 

sup  - 

\\wm\\x 


(23) 


is  achieved  for  wm  =  llvfjij/]  and  thus 


Pm  -  liniv„[(A]||,Y- 


From  relations  (23)  and  (24),  we  infer  that  Hwm  [*A]  is  such  that 


(24) 


\\Jm  [nty„[iA]]  lU  _  J_  <  IIJm  M  Ik 
linivM['A]IU  Pm  ip^  II^IU 

Relations  (22)  and  (25)  yield  the  final  equality  —  =  Km- 

Pm 


(25) 


□ 


Remark  2.5.  The  link  between  the  Lebesgue  constant  Km  and  the  inf-sup  quantity  Pm  introduced  in  Theorem  2.4 
shows  that  Km  depends  on  the  dictionary  of  linear  functionals  E  and  also  on  the  interpolating  space  Xm-  Although 
no  theoretical  analysis  of  the  impact  of  these  elements  has  been  possible  so  far,  we  present  in  Section  4  a  numerical 
study  about  the  influence  of  the  dictionary  E  in  Km- 


Remark  2.6.  Note  that,  since  Theorem  2.4  holds  only  in  Hilbert  spaces,  formula  (16)  does  not  apply  to  the  Lebesgue 
constant  of  the  classical  EIM  given  that  it  is  defined  in  the  L”(Q)  norm.  The  Hilbertian  framework  allows  nevertheless 
to  consider  Dirac  masses  as  linear  functionals  like  in  EIM  if  we  place  ourselves,  e.g.,  in 


2.3.  The  Greedy  algorithm  aims  at  optimizing  the  Lebesgue  constant 

If  we  look  in  detail  at  the  steps  followed  by  the  Greedy  algorithm,  once  Xm-\  and  Wm-\  have  been  derived,  the 
construction  of  Xm  and  Wm  starts  by  adding  an  element  ip  to  Xm-\.  In  the  Greedy  process,  this  is  done  following 
formula  (4),  but  let  us  analyze  what  happens  when  we  add  any  <pm  6  F-  The  first  consequence  of  its  addition  is  that 
the  resulting  inf-sup  constant  becomes  zero: 


(y,  w)x 

(pespanriM_i,i^M)  wewh-i  ll',£’llA'l|vt'IU 


inf 


sup 


=  0. 


(26) 


Indeed,  the  addition  of  pm  to  the  interpolating  basis  functions  has  the  consequence  of  adding  the  element  ipM  - 
<Tm  -  Jm-\Vpm}  that,  by  definition,  satisfies  (ipM,  'iv)x  -  0,  Vw  e  Wm-\-  We  thus  need  to  add  an  element  to  Wm-\  in 
order  to  stabilize  the  inf-sup  condition. 
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Let  us  denote  by  W  the  set  of  Riesz  representations  in  X  of  the  elements  of  our  dictionary  S.  Since 

.  .  (v,  w)x 

inf  sup  - 

(^espan{XM_,,(£.M)  „eWM-\  ll'/’IUI|vv|U 

is  reached  by  the  aim  is  to  add  an  element  wm  of  W  that  maximizes 


(ifiM,  W)x 
v/eW  ||w||x 


(27) 


Since  the  elements  of  the  dictionary  are  of  norm  1  (see  property  PI  above),  this  corresponds  exactly  to  one  of  the  steps 
performed  by  the  Greedy  algorithm  (see  equation  (5)).  Furthermore,  from  the  unisolvence  property  of  our  dictionary, 
the  application 


(f  i->  max  {ip,  w)x 

weW 

defines  a  norm  in  X.  Then,  formula  (27)  reads: 

iiPM,w)x_ 

max — - — -  -  \\ipM  -  SFm-iWmWU- 

weW  ||w||x 

It  is  thus  clear  that  the  choice  of  ipM  that  maximizes  the  value  of  /I is  the  one  that  maximizes  Pm  -  in  the 

ll-llt  norm.  However,  since  in  practice  we  do  not  have  access  to  the  entire  knowledge  of  this  norm,  ||.||,  is  replaced  by 
the  ambient  norm  ||.||;^: 

ipM  =  argmaxll^  -  J'm-iMII.  ~  argmax||^  -  J'm-iMIU,  (28) 

ipeF  ipeF 

which  is  exactly  what  the  Greedy  algorithm  does  (see  (4)).  Hence,  as  a  conclusion,  with  the  practical  tools  that  can  be 
implemented,  the  choice  of  tpM  aims  at  minimizing  the  Lebesgue  constant  with  the  approximation  explained  in  (28). 


3.  Practical  implementation  of  the  Greedy  algorithm  and  the  Lehesgue  constant 

In  the  present  section,  we  discuss  some  practical  issues  regarding  the  implementation  of  the  Greedy  algorithm  and 
the  Lebesgue  constant  Am- 

Since  the  cardinality  of  F  is  usually  infinite,  the  practical  implementation  of  the  Greedy  algorithm  is  carried  out  in 
a  large  enough  sample  subset  Sf  of  finite  cardinality  #Sf  much  larger  than  the  dimension  of  the  discrete  spaces  Xm 
and  Wm  we  plan  to  use.  For  example,  if  F  =  {u{pi,  ■),pi  e  £)),  we  choose  Sf  =  {m(jU,  .),/U  6  c  £))  and  consists 
of  #Sf  parameter  sample  points  pi.  We  assume  that  this  sample  subset  is  representative  enough  of  the  entire  set  F  in 
the  sense  that 

supt  inf  ||x-y||x[ 

xeF  kvespan{.Sf)  J 

is  much  smaller  than  the  accuracy  we  envision  through  the  interpolation  process.  This  assumption  is  valid  for  small 
dimension  of  F,  or,  more  precisely,  for  small  dimension  of  the  parameter  set  D.  In  case  it  cannot  be  implemented 
directly,  we  can  follow  two  strategies  that  have  been  introduced  on  greedy  approaches  for  reduced  basis  approxima¬ 
tions  either  based  on  (parameter)  domain  decomposition  like  in  [15]  or  [16]  based  on  an  adaptive  construction  of  the 
sample  subset,  starting  from  a  very  coarse  definition  as  in  [17].  These  approaches  have  not  been  implemented  here 
but  we  do  not  foresee  any  difficulty  in  adopting  them  to  the  GEIM  framework. 

The  following  lemma  shows  that  the  generalized  interpolant  can  be  recursively  computed. 

Lemma  3.1.  For  any  function  f  €  X,  we  have  the  following  recursion  for  M  >  I 

I J'm[/]  =  J'm-1  [/]  +  CTMif  -  jM-df])qM 

\jo[f]  =0  ^  ^ 

and  the  generalized  interpolant  off  can  be  recursively  computed. 


10 


Proof.  Using  the  fact  that  the  spaces  Xm  are  hierarchically  defined,  both  sides  of  (29)  belong  to  Xm-  Using  the  fact 
that  CTiiqM)  -  0  for  i  <  M  and  the  definition  of  and  we  infer  that 

CT;  (Jm[/])  =  o-i  (Jm-1  [/]  +  o-M  (/  -  Jm-\  [/])  qm) , 

Finally,  it  is  clear  that  the  right  and  left  hand  sides  have  the  same  image  trough  ctm-  The  equality  holds  by  uniqueness 
of  the  generalized  interpolation  procedure. 

□ 


Remark  3.2.  This  result  also  holds  for  the  classical  EIM  case. 

The  greedy  algorithm  is  in  practice  a  very  time-consuming  task  whose  computing  time  could  significantly  be 
reduced  by  the  use  of  parallel  architectures  and  the  use  of  formula  (29)  as  is  outlined  in  Algorithm  1 . 


Algorithm  1  Practical  implementation  of  the  Greedy  procedure 

1 

Input:  E,  Sf  =  [fk  6  T)*ff ,  e,„i,  M„ax,  M  =  0 

2 

Assign  a  set  of  functions  ...,  fkp^^op]  to  each  processor  p. 

3 

repeat 

4 

M  M+  1 

5 

>  parallel 

6 

for  k  =  {kp  start-:  •  •  ■ » kp  stop]  do 

7 

f  =  fk 

8 

Compute  and  store  a-^f/  -  STmif))- 

9 

Assemble  J7jm+i(/)  following  formula  (29) 

10 

Compute  Sm+1  =  \\f  -  Jm+ADWx 

11 

if  >■  Sp  Ytiax  then 

12 

kp,msL\  ~  ^  £p,rnax  ~  ^M-i-1 

13 

end  if 

14 

end  for 

>  end  parallel 

15 

Gather  max)}  find  (fTjjjax)  ^max)  “  arg  HiaX  (£^p,maxi  ^/7,max)- 

^  pE\l,...,Nproc\ 

16 

^M+l  fkmax  ^^./^max^ 

17 

^/>,max  “  0 

>  parallel 

18 

thr  j  —  {jp^startf  •  •  ■  >  jp^stop] 

19 

cr  =  CTJ 

20 

Compute  Sm+1  =  |cr(rM+i)l 

21 

if  £m+i  ^  ^/>,max  then 

22 

jp.mnK  ~  j  Snd  £y7,max  ~  ^M-i-1 

23 

end  if 

24 

end  for 

>  end  parallel 

25 

26 

(  I  proc 

Gather  1  (fip^tnax?  7/7, max)|  ,  ^d  find  (fijnaxi  y'max)  ~  arg  maX  maxi  max)- 

pe[U...,Nproc\ 

Compute  and  store  qm+i  =  - -7— — r  ■ 

27 

Store  o-M+i  =  cTj^. 

28 

Compute  and  store  wm+i  (Riesz  representation  of  ctm+i)- 

29 

until  Smax  <  e,„i  01 M  >  M^ax 

30 

Output:  {a-,,...,o-M+il,  Wm+i  =  span{wi, . . . ,  wm+iI,  Xm+i  =  span{^i, . . . ,  ?m+iI- 

Once  Xm  and  Wm  have  been  constructed  thanks  to  Algorithm  1,  the  Lebesgue  constant  can  be  computed  by  the 
resolution  of  an  eigenvalue  problem  as  is  explained  in 

Lemma  3.3.  If{qi, . . . ,  §m)  ond  {vvi, . . . ,  wm}  are  an  orthonormal  basis  ofXu  and  Wm  respectively,  then 

Pm  =  1/Am  =  (30) 

where  A  is  the  M  X  M  matrix  whose  entries  are  Aij  —  (vv,-,  qpx  and  A^niaiA^A)  denotes  the  minimum  eigenvalue  of  the 
positive  definite  matrix  A^ A. 
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Proof.  Since 


„  .  .  (x,w)x 

Pm  =  ini  sup - 


.  {Ax,  w)2 

=  inr  sup - 

u;eR«ll-*ll2l|w||2 


.  .  IIA^Ib 

inf - , 

A:eB«  ||x||2 


the  result  easily  follows  because 


IIAxll^ 

l|x||2 


is  the  Rayleigh  quotient  of  A^A  whose  inhmum  is  achieved  by  A^in{A^A). 


□ 


Remark  3.4.  Note  that  Pm  corresponds  to  the  minimum  singular  value  of  the  matrix  A,  which  is  a  matrix  of  small 
size  M  X  M.  Its  computation  can  be  easily  performed  by,  e.g.,  the  inverse  power  method. 


4.  A  numerical  study  about  the  impact  of  the  dictionary  S  of  linear  functionals  in  the  Lebesgue  constant 

As  outlined  in  Remark  2.5,  the  explicit  expression  of  the  Lebesgue  constant  presented  in  formula  (16)  shows  that 
Am  is  intimately  linked  to  the  dictionary  of  linear  functionals  S  that  is  used  in  the  Greedy  algorithm  to  build  the 
interpolation  process.  With  the  exception  of  the  trivial  case  considered  in  Corollary  2.2,  no  theoretical  analysis  of 
the  impact  of  2  on  the  behavior  of  the  Lebesgue  constant  has  been  possible  so  far.  For  this  reason,  we  present  here 
some  numerical  results  on  this  issue  as  a  first  illustration  of  this  connection.  The  same  computations  will  also  let  us 
numerically  validate  the  formula  (16)  for  Am,  whose  original  definition  is  given  by  (7). 

We  place  ourselves  in  Q  =  [0,1]  and  consider  the  numerical  approximation  in  L^(Q)  or  H^{Q)  of  the  following 
compact  set: 


F  =  {/(.,P1,P2)I(//1,JU2)6  [0.01,24.9]X[0,15]),  (31) 

where 

f{x,pi,p2)  ^  ^  VxeG. 

yl  +  (25  +  Pi  cos{p2x))x^ 

We  recall  that  -  {/  |  ||/||£2(Q)  <  co),  where  the  norm  ||  ■  ||L2(n)  is  induced  by  the  inner  product  (w,  v)£2(Q)  = 

J  w(x)v(x)dx.  Also,  //*(f2)  =  {/ 1  ||/||/£i(n)  <  oo),  where  the  norm  ||  ■  llnqn)  is  induced  by  the  inner  product  {w,  v)/£i(n)  = 
n 

f  w(x)v(x)dx  +  J  Vw(x).Vv(x)dx. 
n  n 

Any  f  &  F  will  be  approximated  by  its  generalized  interpolant  at  dimension  M.  For  this  purpose,  the  practical 
construction  of  the  interpolating  space  Xm  and  the  selection  of  the  linear  functionals  is  done  through  the  Greedy 
algorithm  described  in  Section  3.  The  following  dictionary  of  linear  functionals  has  been  employed: 


S  =  (o-^  6  Z{X),  ke{l,.. . , 


(32) 


where  N sensor  -  150,  and 


The  function  q  ^  reads: 


where 


o-kitp)  =  J"  Ck^s{x)ip{x)dx, 

.xeSl 


Ck,s{x)  = 


mk.six) 


INA:,i(OllL>(n) 

mk,six) 


SipeX. 

VxeQ, 

VxeQ 


(33) 


and  Xk  e  G.  We  will  explore  the  variation  of  the  coefficient  s  e  IR+  in  order  to  understand  the  influence  of  the 
dictionary  E  on  A^. 
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4.1.  Validation  of  the  inf-sup  formula 

We  will  first  start  by  fixing  s  to  a  value  of  0.005  and  by  numerically  validating  formula  (16)  of  the  Lebesgue 
constant  by  comparing  it  to  the  value  given  by  the  original  formula  (7). 

Regarding  the  computation  of  (16),  the  quantity  /Sm  has  been  derived  using  formula  (30)  of  Lemma  3.3.  It  suffices 
to  evaluate  the  scalar  products  of  the  matrix  A  defined  in  that  lemma  and  obtain  the  minimum  eigenvalue  of  A^A. 
For  the  practical  computations,  a  Pi  finite  element  approximation  of  the  functions  q,  and  w,  has  been  used  in  order 
to  simplify  the  scalar  product  evaluation  in  the  L?  and  spaces.  For  the  same  reason  and  as  a  matter  of  global 
coherence,  the  Lebesgue  constant 

IIJ'mMIU 

(peX  ll^lw 

is  also  approximated  in  the  same  Pi  finite  element  approximation  of  the  elements  of  ?(.  This  approach  leads  to  the 
computation  of  a  discrete  Raleigh  quotient,  whose  derivation  is  explained  in  detail  in  appendix  B. 

The  results  of  the  computation  are  given  in  Figure  2  and  show  an  excellent  agreement  between  both  values  in 
and  //' .  The  same  agreement  holds  for  any  value  of  the  parameter  s  of  the  linear  functionals,  but,  as  will  be  presented 
in  the  next  section,  the  behavior  of  Am  varies  depending  on  this  parameter. 


10 


10" 


10 


107 


10' 


I  ‘  1 

o  A. .  with  sup  formula 
M  ^ 

o  A. .  with  sup  formula 

1  M  ^  1 

^  2 
10 

< 

10 


10  20  30  40 
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50 


10' 


%8 


0 


(a)A  =  L2([0;l])  (b)  A  = //‘([O;  1]) 

Figure  2:  Numerical  validation  of  the  inf-sup  formula:  comparison  between  formulae  (7)  and  (16). 


10  20  30  40  50 

Dimension  M 


In  the  particular  case  presented  here,  the  behavior  of  the  Lebesgue  constant  does  not  significantly  change  if  we 
place  ourselves  in  L?  or  in  //'  and  Am  remains  constant  (the  degradation  in  the  behavior  for  M  >  44  is  due  to 
numerical  round-off  errors). 


4.2.  Impact  of  the  dictionary  of  linear  functionals 

We  now  study  the  impact  of  s  on  the  evolution  of  the  Lebesgue  constant  through  our  example  in  one  dimension. 
For  this  purpose,  we  present  in  Figures  3a  and  3b  the  behavior  in  and  in  of  Am  for  different  values  of  s. 

To  begin  with,  we  will  focus  on  the  behavior  for  sufficiently  large  values  of  s  and  analyze  the  range  s  >  5.10*^. 
It  can  be  observed  that,  as  s  increases,  the  Lebesgue  constant  progressively  degrades  in  both  norms.  The  sequence 
(Am)  starts  to  diverge  at  dimensions  that  are  lower  and  lower  as  s  increases  (compare,  e.g.,  the  behaviors  between  the 
case  s  -  2.10“^  and  s  -  4.10“^).  An  intuitive  manner  to  interpret  this  observation  is  as  follows;  the  dictionary  under 
consideration  in  this  example  (see  formula  (32))  consists  on  local  averages  operations  whose  ’’range”  is  controlled  by 
s.  As  s  increases,  the  range  increases  and  a  limit  will  be  reached  in  which  the  addition  of  more  linear  functionals  will 
result  in  a  redundant  addition  of  information  because  of  an  overlap  of  the  domains  where  the  local  averages  are  acting. 
As  a  result,  the  larger  s,  the  sooner  this  redundancy  will  appear  and  the  more  unstable  the  process. 
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Dimension  M 


Dimension  M 


(a)L"  norm. 


(b)//‘  norm. 


Figure  3:  Impact  on  Am  of  the  parameter  s  of  the  linear  functionals. 


It  is  also  important  to  understand  the  behavior  when  the  parameter  s  tends  to  zero.  In  this  case,  the  linear  func¬ 
tionals  tend  to  Dirac  masses,  that,  in  ID,  are  elements  of  //  '  but  not  of  L^.  Hence,  in  the  limit  i  =  0,  the  dehnition  of 
the  space  Wm  will  be  possible  in  //'  but  not  in  because  the  problem; 


Find  w,  e  X  such  that; 

o-iiip)  =  (w;,  (p)x  =  'iipeX 


(34) 


is  well-dehned  in  i/'  but  not  in  L^.  This  observation  helps  to  understand  first  of  all  why  Ai  remains  roughly  constant 
in  //'  as  s  decreases  whereas  it  behaves  as  in  the  l}  norm  (see  Figure  4).  Indeed,  in  the  //'  case,  we  have  the 
inequality 


11.71  Mll/rqn)  ,  .  ..II^iIIhho)  ^  ll?illHi(n) 

— pn - =  ki(<^)|-— r - <  IMlL”(n)7p - , 

ll'/'ll/rqn)  IIv^IIhico)  II^^IIhh^) 


e 


which  is  bounded  for  any  i  e  IR+.  However,  in  the  case  of  it  can  be  infeiTed  that 


\VI\[^]\\l'^(SI)  .  ,  ,,  II^iIIl^IO) 

— fTi - "  —  - 


INi..9llL^(n) 

INi,.sllz,>(n) 


ll?illL2(n), 


€  h\q.) 


X  —  X] 

where  we  have  applied  the  Cauchy-Schwarz  inequality  to  \a-\{ip)\.  A  simple  change  of  variable  u  -  -  in  the 

s 

\\mi  5lh2(o^ 

evaluation  of - ^ - -  leads  to  the  bound 

INl,.sllL>(ti) 


WJimrnn) 


V(^  e  L\n), 


(35) 


where 

f  e~‘‘^du 

C  -  — _ . 

J  g-u-jldii 

n 

In  Figure  4,  note  that  for  values  s  <  10“^,  the  behavior  of  Ai  no  longer  follows  but  this  is  due  to  computer 
limitations.  Indeed,  the  computations  have  been  carried  out  with  a  maximum  number  of  10"^  degrees  of  freedom  in 
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the  Pi  approximation  because  of  memory  storage  issues.  As  a  result,  for  s  <  10  we  no  longer  capture  enough 
information  with  this  hnite  element  precision. 


Figure  4:  Behavior  of  Ai  as  a  function  of  i  (//'  and  norms).  Remark:  the  scale  of  the  figure  is  log-log. 


As  a  consequence  of  the  diverging  behavior  of  Ai  in  as  the  parameter  s  decreases,  it  is  reasonable  to  expect 
that  the  sequence  (Am)  quickly  diverges  as  s  ^  0  in  but  that  it  remains  bounded  in  //'.  This  behavior  is  indeed 
illustrated  in  Figures  3a  and  3b  through  the  example  of  i  =  10  in  which  it  is  possible  to  observe  the  phenomenon. 

5.  Application  of  GEIM  to  the  real-time  monitoring  of  a  physical  experiment 

The  main  purpose  of  this  section  is  to  illustrate  that  GEIM  can  be  used  as  a  tool  for  the  real-time  monitoring 
of  a  physical  or  industrial  process.  The  rationale  is  to  provide  an  online  accurate  approximation  of  the  phenomenon 
under  consideration  thanks  to  the  computation  of  a  generalized  interpolating  function  that  will  be  derived  on  the  fly  by 
combining  in  an  appropriate  manner  measurements  from  the  experiment  with  mathematical  models  (a  parameter  de¬ 
pendent  PDE).  Such  a  tool  could  be  employed  for  keeping  track  of  the  process  in  the  whole  domain  of  the  experiment 
and  not  only  at  locations  where  the  sensors  are  placed. 

We  will  also  show  that  the  proposed  method  can  be  helpful  in  design  of  experiments:  the  method  minimizes  the 
number  of  sensors  required  to  reconstruct  the  field  variable  and  informs  the  optimal  sensor  types  and  placement. 

5.7.  The  general  method 

Let  us  assume  that  we  want  to  monitor  in  real  time  a  field  Mtme  appearing  as  an  input  for  some  quantities  of 
interest  in  a  given  experiment  that  involves  sensor  measurements.  We  assume  that  the  conditions  of  the  experiment 
are  described  by  a  vector  of  parameters  //tme  6  E,  where  £  is  a  compact  set  of  and  that  Mu-ue  is  the  solution  of  a 
parameter  dependent  PDE 

Df,u  ^  g^  lie  E  (36) 

when  It  =  //true  (in  other  words  Mtiue  =  Ths  vector  //true  will  be  unknown  in  general  so  the  computation  of 

Mtrue  cannot  be  done  by  traditional  discretization  techniques  like  hnite  elements.  Besides,  even  if  //tme  was  known, 
its  computation  could  not  be  performed  in  real-time  with  classical  techniques.  For  all  these  reasons,  we  propose  to 
compute  the  generalized  interpolant  as  an  approximation  of  i/tme- 

Such  an  approximation  requires  that  the  set  of  solutions  {m^,  V//  6  E]  is  included  in  some  compact  set  7^  of  A  that 
is  of  small  Kolmogorov  n-width  in  X  ([12]).  A  dictionary  S  c  JL{X)  is  also  required.  Each  element  cr  e  -L{X)  will 
mathematically  model  a  sensor  in  the  experiment.  The  location  and  type  of  the  sensors  might  be  hxed  a  priori  (by 
the  experts  running  the  experiment),  but  it  is  also  possible  to  consider  that  2  is  composed  by  elements  that  represent 
potential  locations  and  potential  types  of  sensors. 
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Since  we  need  to  define  the  generalized  interpolating  spaces  Xm  -  span{^i, . . .  ,qM}  together  with  the  suitable 
interpolating  linear  functionals  {cri, . . . ,  cr^),  a  greedy  algorithm  has  to  be  performed  beforehand  and  therefore  the 
computation  of  is  divided  into  two  steps: 

•  In  an  offline  phase  (i.e.  before  the  experiment  takes  place): 

-  We  define  a  finite  subset  Sf  -  {uip,  .),p  e  (Z  E}  c  F  and  solve  (36)  for  each  element  of  Sf  with  an 
accurate  enough  discretization  strategy.  This  can  be  done  with  traditional  approximation  tools  like,  e.g., 
finite  elements  or  a  reduced  basis  strategy. 

-  Following  the  steps  of  Algorithm  1,  a  greedy  algorithm  over  the  set  Sf  performed  to  build  an  M- 
dimensional  reduced  basis  Xm  =  span{q'j  e  F,j  6  [1,M])  together  with  the  suitable  linear  functionals 
{(Ti, . . . ,  ctm}  coming  from  the  dictionary  2.  The  selection  of  the  linear  functionals  means  that,  among  all 
the  sensors  in  the  experiment  that  constitute  our  dictionary  S,  we  select  the  M  most  suitable  according 
to  the  greedy  criterion.  These  M  linear  functionals/sensors  will  be  the  only  ones  required  to  reconstruct 
the  field  so  the  present  methodology  could  be  viewed  as  a  strategy  to  minimize  the  number  of  sensors 
involved  in  the  experiment  as  well  as  to  find  their  optimal  type  and  placement. 

•  In  an  online  phase  (i.e.  when  the  experiment  is  running),  we  collect  in  real  time  the  measurements 

from  the  M  selected  sensors  (that  are  placed  in  the  experiment).  The  generalized  interpolant  can  then 

be  computed  following  formula  (3).  It  has  been  observed  so  far  (see  the  numerical  example  below  and  [1]) 
that  the  interpolation  error  decreases  very  quickly  as  the  dimension  M  increases  and  therefore  relatively  small 
values  of  M  are  required  to  reach  a  good  accuracy  in  the  approximation  of  by  Thanks  to  this, 

the  computation  of  can  be  performed  in  real-time  (or  almost). 

Remark  5.1.  Note  that  our  strategy  supposes  that  the  physical  experiment  u,rue  is  perfectly  described  by  the  solution 
u^  of  (36)  for  p  =  Ptme-  This  is  a  very  strong  hypothesis  because  the  model  might  not  perfectly  describe  the  experiment 
under  consideration.  Besides,  it  is  here  assumed  that  there  is  no  noise  in  the  measurements,  which  is  also  a  strong 
assumption.  In  [1  ],  some  preliminary  analysis  has  been  presented  to  take  into  account  the  presence  of  noise  in  the 
measurements.  Regarding  the  model  bias,  in  the  recent  works  of  [18,  19],  the  authors  are  able  to  take  it  into  account 
under  several  hypothesis  in  the  so  called  "Parametrized-Background  Data-Weak  Formulation”  for  variational  data 
assimilation.  In  fact,  GEIM  is  a  particular  instance  of  this  method  for  the  case  (with  the  notations  of  [19])  N  —  M 
and  this  latter  choice  is  appropriate  for  situations  in  which  the  bias  is  small. 

Remark  5.2.  In  the  strategy  proposed  in  this  section,  sensor  measurements  are  incorporated  in  the  interpolation 
procedure  through  the  space  Wm  (which  is  spanned  by  the  Riesz  representations  of  the  linear  functionals  of  the 
sensors).  In  the  reference  [20],  one  can  find  an  early  work  in  oceanography  in  which  data  assimilation  is  also 
incorporated  through  the  construction  of  the  space  Wm-  However,  in  the  case  of  [20],  no  a  priori  error  analysis  was 
provided  in  the  computational  procedure  that  was  proposed. 


5.2.  A  numerical  application  involving  the  Stokes  equation 

We  are  going  to  illustrate  the  procedure  in  the  case  where  the  experiment  corresponds  to  a  lid-driven  cavity 
problem  that  takes  place  in  the  spatial  domain  Q.  -  [0;  1]  x  [0;  1]  c  We  consider  two  parameters  p  -  (qi\,p2)  6 
[1;  8]  X  [1;  8]  such  that,  for  a  given  p,  the  parametrized  PDF  reads 


Find  the  solution  {Uij,pf)  6  (//’(Q))  x  Lq(Q.)  of  : 
-AUfj  -H  grad(p^)  =  f^,  a.e.  in  Q 


div(uM)  =  0,  a.e.  in  Q 


x(l  -  x) 

0 


,  a.e.  on  Fi 


V  ”  / 

=  0,  a.e.  on  80.  \  F i 


(37) 
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where  the  forcing  term  = 


100sin(//iny) 

-lOOsin 


and  Fi  =  {x  e  [0;  1],  y  -  1).  The  space  Li{Q.)  corresponds  to  the 


functions  in  L'^(Q)  with  zero  mean  over  Q,  namely 


l2(Q)  =  |/ 6  J^/  =  o|. 

Two  examples  of  solutions  are  provided  on  Figures  5  and  6  where  finite  elements  ^  Pl-bubble/Pl  have  been  used  to 
approximate  the  velocity/pressure  fields  respectively  (this  approximation  is  the  one  that  will  be  used  to  compute  the 
off-line  snapshots  that  form  Sf). 

■  0.456 

-0.053 
1-0.561 

Figure  5:  From  left  to  right:  pressure,  horizontal  and  vertical  velocity  solutions  for  the  parameter  ^  =  (5;  1) . 


Figure  6:  From  left  to  right:  pressure,  horizontal  and  vertical  velocity  solutions  for  the  parameter  fi  =  (8;  5) . 


We  assume  that 

•  the  set  of  solutions  {(u,p)(//),  Vyu)  c  F  and  F  is  of  small  Kolmogorov  n-width  in  X  Lq(0).  This 

assumption  is  made  a  priori  and  will  be  verified  a  posteriori  in  a  convergence  study  of  the  interpolation  errors. 

•  we  have  velocity  and  pressure  sensors  at  our  disposal  which  mathematically  means  that  we  have: 

-  a  dictionary  for  the  velocity:  2“  =  {cr“)  c 

-  a  dictionary  for  the  pressure:  =  {erf)  c 

In  our  numerical  example,  the  linear  functionals  that  have  been  used  consist  of  local  averages  of  the  same  form  as 
(32)  and  (33)  but  adapted  to  the  2D  case.  The  parameter  s  has  been  fixed  to  i  =  10~^  and  we  will  have  N sensor  -  100 
sensors  for  the  pressure  and  other  N sensor  -  100  sensors  for  the  velocity.  The  centers  of  these  local  averages  are 
located  on  a  10  x  10  equispaced  grid  over  D. 

Given  an  experiment  corresponding  to  the  vector  of  parameters  jUexp,  we  are  going  to  —  quickly  and  accurately — 
approximate  in  G  the  vectorial  field  (u,  p)(pexp)  by  its  generalized  interpolant  J'm  [(u,  p)(juexp)]  using  the  only  knowl¬ 
edge  of  measurements  from  sensors.  Because  we  are  facing  here  the  reconstruction  of  a  vectorial  field,  several 
potential  input  from  (u,  p){p)  can  be  proposed.  In  the  present  paper,  three  different  classes  will  be  considered.  They 
will  all  fulfill  the  divergence-free  condition  for  the  velocity  interpolant  div  {J'm  [u(A()])  =  0. 

Reconstruction  1:  Independent  treatment  of  u(p)  and  p{p). 

The  first  possibility  is  to  treat  (u,  p){p)  not  as  a  vectorial  field  but  as  two  independent  fields  u(p)  and  p{p)  to  interpolate 


*Our  computations  have  been  carried  out  with  Freefem++[21],  [22], 
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independently  with  velocity  measurements  for  u(p)  and  pressure  measurements  for  pijj).  In  other  words,  the  gener¬ 
alized  interpolant  is  defined  in  this  case  as  [(Ujf’)(A')]  =  {<!Tm  [p(/t)])-  This  requires  the  offline 

execution  of  two  greedy  algorithms;  one  for  the  velocity  and  another  for  the  pressure.  Each  one  respectively  provides 

•  a  velocity  basis  {u(ju,))|^j  and  a  set  of  Mu  velocity  sensors  {cr“)|^j  chosen  among  the  dictionary  E“.  The  inter- 

Mu 

polant  for  the  velocity  will  be  [u(ju)]  =  2  aiUijUj)  where  the  o',  are  given  by  the  interpolating  conditions 

"  !=1 

JuOu)])  =  cr“(uOu)),  Vf  e  {1, . . . ,  Mu). 

•  a  pressure  basis  ^  set  of  pressure  sensors  chosen  among  The  interpolant  for  the 

Mp 

pressure  will  be  -  2  7jP(Pj)  where  the  yj  are  given  by  the  interpolating  conditions  [p(p)])  - 

J\  p  } 

cr'’{p(p)),  V;e{l,...,Mp). 

Note  that  in  this  approximation,  the  construction  of  p)ip)\  involves  Mp  pressure  sensors  and  Mu  veloc¬ 

ity  sensors,  i.e.  Mp  +  Mu  coefficients.  In  Figure  7,  we  have  represented  the  locations  of  the  sensors  in  the  order  given 
by  the  greedy  algorithm. 


1 

’ 

’ 

1 

0.9 

□  02 

□  10 

□  25 

□  04 

□  06 

□  08 

□  05 

0.9 

□  10 

□  23 

□  20 

□  17 

□  15  - 

0.8 

n  14 

□  19 

□  20 

□  30 

0.8 

D  14 

□  26 

□  24  - 

0.7 

□  17 

□  27 

□  29 

0.7 

□  18 

□  12 

□  29 

0.6 

□  07 

□  18 

0.6 

□  28 

□  13 

□  04 

• 

0.5 

□  21 

0.5 

□  11 

□  02 

□  05 

□  30  ^ 

0.4 

□  22 

□  03 

0.4 

□  03 

□  06  - 

0.3 

□  12 

□  28 

0.3 

□  08 

□  27 

0.2 

□  15 

□  24 

□  16 

0.2 

□  09 

□  07 

□  01 

□  19  ^ 

0.1 

□  23 

□  09 

□  13 

□  26 

□  11 

□  01 

0.1 

B  16 

□  25 

□  22 

□  21  ^ 

0.2 

0.4 

0.6 

0.8 

1  "o 

0.2 

0.4 

0.6 

0.8 

(a)  Pressure  sensors  (b)  Velocity  sensors 

Figure  7:  Locations  of  the  sensors  for  reconstruction  1. 


The  performances  of  the  method  are  plotted  in  Figure  8  where  a  numerical  estimation  of  the  behavior  of  the 
interpolating  errors  for  the  reconstruction  of  u  and  p  have  been  represented.  These  values  have  been  obtained  by  the 
interpolation  of  196  configurations  coming  from  different  parameter  values  pi  following  formula; 


\\p(pd  -  XplP(^^i)]\\LHn) 

fed,..., 196)  WpipdWLHn) 

[uOu/)]  IlHffnp 

/e|l....,196)  l|uOu/)||Hi(np 


(38) 


In  this  figure,  we  can  observe  the  convergence  of  the  interpolation  errors  for  both  the  velocity  and  pressure  fields. 
After  a  preasymptotic  stage  for  interpolating  spaces  of  small  dimension,  an  exponential  convergence  of  the  error  is 
observed.  After  about  dimension  M  =  25,  the  error  stagnates  due  to  the  fact  that  we  have  reached  the  finite  element 
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accuracy  used  for  the  computation  of  the  offline  snapshots.  The  computation  of  the  Lebesgue  constants 


M„ 


lIJ'f 


sup  - 

peLHSl) 


M. 


(p)\y 


(it) 


IlfllL^cn) 


K,  =  sup 

ueH'iClf 


lIJ'M.WIlHhny 

l|u||Hi(n)^ 


(39) 


has  also  been  performed  following  formula  (16).  Its  behavior  seems  linear  with  the  dimension  of  interpolation  and  is 
therefore  far  from  the  crude  theoretical  upper  bound  given  in  formula  (8).  From  the  results  presented  in  Section  4,  an 
idea  to  improve  the  behavior  of  the  Lebesgue  constant  could  be  to  consider  a  smaller  value  for  s.  However,  this  does 
not  seem  necessary  here  given  the  good  convergence  properties  observed.  Furthermore,  in  a  real  case,  the  parameter 
s  of  the  linear  functionals  is  difficult  to  change  as  it  is  fixed  by  the  filter  characteristics  of  the  sensors  involved  in  the 
experiment. 

Remark  5.3.  In  the  case  of  a  quickly  diverging  Lebesgue  constant,  the  convergence  properties  of  the  method  could  be 
severely  degraded  (see  Theorem  1.4).  Ai  Section  4  illustrates,  one  of  the  reasons  for  the  divergence  could  be  related 
to  a  lack  filter-width  of  the  sensors  (parameter  s  too  large)  and,  in  this  case,  it  would  be  necessary  to  restart  the 
procedure  with  sensors  with  a  smaller  filter-width. 


Figure  8:  Reconstruction  1:  A  numerical  estimation  of  the  behavior  of  the  interpolation  error  (left)  and  the  Lebesgue  constant  (right)  as  a  function 
of  the  dimension  of  the  interpolating  spaces  Xm 


Reconstructions  2  and  3:  Vectorial  treatment  for  u(p)  and  p(p). 

An  alternative  to  the  first  reconstruction  is  to  consider  (u,  p)ip)  as  a  vectorial  field  and  define  its  generalized  interpolant 
as  ffu  [(u,p)(/t)]  :=  Yjfii  Ji  (u>f)  (F/)’  where  now  only  M  coefficients  7,  are  involved.  The  joint  basis  {(u,p)  (/t,))^[ 
is  provided  by  a  greedy  algorithm  in  the  offline  stage  together  with  a  set  of  M  linear  functionals  Each 

of  these  linear  functionals  involve  pressure  and  velocity  measurements  at  a  given  spatial  location  and  are  defined  as 
:=  cr"(u)  +  cr^(p).  The  interpolating  conditions  for  the  inference  of  the  coefficients  y,  are  now  the  following; 

M 

p)(p))  =  [(u,  p)(p)]  )  =  ^  (u,  p)  (pj)),  Vi  e  { 1 , . . . ,  M),  (40) 

>=i 

Notice  that  this  definition  of  the  linear  functionals  can  involve  both  velocity  and  pressure  measurements  or 

can  take  into  account  velocity  or  pressure  measurements  only  by  setting  cr“  =  0  or  cr'’  =  0.  We  have  explored  this 
flexibility  in  the  following  two  reconstructions; 

•  the  interpolation  of  the  pressure  and  velocity  fields  with  pressure  and  velocity  measurements;  ;=  cr“(u)  + 
cr^j(p)  (reconstruction  2). 
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•  the  interpolation  of  the  pressure  and  velocity  fields  with  pressure  measurements  only;  ;=  cr^ip).  In  other 
words,  we  are  here  studying  if  a  velocity  field  can  efficiently  be  reconstructed  with  the  only  knowledge  of 
pressure  measurements  (reconstruction  3). 


The  sensor  locations  provided  by  the  greedy  algorithm  are  shown  in  Figure  9  and  the  results  are  summarized  in 
Figures  10  where  an  estimation  of  the  interpolation  error  is  plotted  according  to  formula  (41). 
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max 


||(u,/7)0u,)  -  [(u,  j7)(/r,)]  ||tfi(n)2xL^(n) 

ll(u,p)(/r,)ll/ri(n)2xL2(n) 


(41) 
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Figure  9:  Locations  of  the  sensors  for  reconstructions  2  and  3. 
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Figure  10:  Reconstructions  2  and  3:  A  numerical  estimation  of  the  behavior  of  the  interpolation  error  (left)  and  the  Lebesgue  constant  (right)  as  a 
function  of  the  dimension  of  the  interpolating  spaces  Xm- 


20 


The  interpolating  error  of  the  two  types  of  reconstructions  presents  a  very  similar  decay  behavior  in  both  cases  and 
the  convergence  is  also  very  similar  to  reconstruction  1 .  The  most  interesting  consequence  of  this  is  that  the  velocity 
can  efficiently  be  reconstructed  with  only  pressure  measurements.  This  result  cannot  probably  be  generalized  to  all 
types  of  situations  but  it  proves  that  in  some  cases  like  the  current  one  there  is  some  redundancy  in  the  data  and  that, 
in  this  precise  problem,  there  is  no  need  in  having  velocity  measurements  in  order  to  obtain  a  good  accuracy  in  the 
approximation  of  the  velocity  field. 

The  Lebesgue  constant 

iu.p)  _ 

~  11/  Ml 

(u,p)eH\n)^xLHQ.)  P)\\H\nfxL'^{n) 

has  also  been  computed  for  reconstructions  2  and  3  as  is  shown  in  Figure  10.  Once  again,  the  behavior  is  linear  which 
is  a  moderate  growth  rate. 

6.  Conclusion  and  perspectives 

After  revisiting  the  foundations  of  GEIM  for  Banach  spaces,  the  present  work  has  focused  on  understanding  the 
stability  of  the  process  and  a  relation  between  Am  and  an  inf-sup  problem  has  been  established  in  the  particular  case 
of  Hilbert  spaces.  An  interpretation  of  the  generalized  interpolant  as  an  oblique  projection  has  also  been  presented 
in  that  case.  The  derived  formula  for  Am  has  also  allowed  us  to  notice  that  the  Greedy  algorithm  optimizes  in  some 
sense  the  Lebesgue  constant. 

A  first  analysis  about  the  impact  of  the  dictionary  of  linear  functionals  L  on  the  Lebesgue  constant  has  also  been 
presented  through  a  numerical  test  case.  Lurthermore,  for  a  given  dictionary  L,  the  Lebesgue  constant  depends  on  the 
norm  of  the  ambient  space  X  (see  formula  (7)).  A  comparison  of  the  behavior  of  (Am)  when  X  -  L}  or  //'  has  been 
provided  in  the  case  of  a  dictionary  composed  of  simple  local  averages. 

Beyond  these  results,  there  are  still  many  challenging  theoretical  open  questions.  Among  the  most  important  we 
mention: 

•  Can  we  obtain  a  general  theory  for  the  impact  of  L  on  the  behavior  of  (Am)?  Can  we  obtain  a  tighter  upper 
bound  than  (8)7 

•  How  can  we  carry  out  the  offline  computation  in  a  reasonable  time  when  the  number  of  parameters  is  large? 

•  How  can  we  treat  the  model  bias  between  Wtrue  the  manifold  of  solutions  of  our  parameter  dependent  PDL? 
The  works  of  [18]  will  probably  be  helpful  to  carry  out  this  task. 

•  How  can  we  treat  noisy  measurements?  One  can  find  some  preliminary  ideas  in  [1]  and  the  works  of  [23]. 

Lurthermore,  the  recent  results  of  [12]  lead  us  to  think  that  it  would  be  interesting  to  explore  non-linear  inputs  of  the 
form 

cr  {t{ip)) , 

where  cr  6  JLiX),  f  :  A  — »  A  is  a  non-linear  mapping  and  ip  is  an  element  of  a  compact  set  of  small  Kolmogorov 
n-width  in  A.  In  an  ongoing  work,  we  are  exploring  this  idea  in  the  case  of  the  Navier-Stokes  equations. 

In  the  second  part  of  the  paper,  we  have  illustrated  one  of  the  most  straightforward  practical  applications  of  GEIM 
that  consists  of  monitoring  in  real-time  a  process.  The  idea  is  that  GEIM  could  reconstruct  in  real-time  physical 
quantities  in  the  whole  domain  of  an  experiment  by  combining  the  real-time  acquisition  of  measurements  from  sensors 
with  mathematical  models  (parameter  dependent  PDEs). 

This  scheme  has  been  applied  to  an  example  dealing  with  a  parametrized  lid-driven  Stokes  equation.  The  example 
shows  a  fast  decrease  in  the  interpolation  etTor,  which  confirms  that  it  is  feasible  to  use  GEIM  to  monitor  experiments 
in  real-time  in  cases  where  dn(F,X)  is  small  enough  (i.e.  when  the  experiment  is  simple  enough).  The  behavior  of 
the  Lebesgue  constant  seems  to  be  linear  and  seems  to  be  in  accordance  with  previous  works  for  the  classical  EIM 
(see  [8]).  The  linear  increase  is  far  from  the  theoretical  exponential  upper  bound  of  (8)  and  suggests  that  the  bound 
might  not  be  optimal  in  sets  F  of  small  Kolmogorov  n-width.  In  the  example,  two  types  of  sensors  have  been  used  (of 
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Experiment 


Figure  11:  A  tool  to  supervise  in  real-time  the  safety  of  an  experiment. 


pressure  and  velocity)  and  the  idea  of  introducing  different  classes  of  sensors  could  be  extended  to  make  more  refined 
distinctions  between  them. 

By  taking  this  method  as  a  starting  point,  GEIM  could  be  used  to  devise  a  more  complete  tool  capable  of  supervis¬ 
ing  the  safety  of  processes  (see  Figure  11).  The  idea  would  be  the  following:  suppose  we  have  2M  measurements  at 
our  disposal.  We  first  invoke  GEIM  to  reconstruct  the  solution  field  using  M  selected  measurements.  We  next  predict, 
based  on  the  GEIM-reconstructed  field,  the  output  associated  with  the  remaining  M  sensors.  We  then  compare  the  M 
predicted  values  with  the  actual  experimental  measurements.  If  the  values  differ  too  much  from  each  other,  then  we 
consider  that  an  abnormal  event  has  occurred  in  the  experiment  and  an  alarm  can  be  sounded  to  signal  the  incident. 

In  addition  to  this  alarm  information,  we  can  seek  to  provide  an  accurate  enough  reconstruction  of  the  solution 
during  the  incident  by  using  the  following  strategy:  through  the  computation  of  an  posteriori  error  estimator  in  the 
regions  where  the  sensor  measurements  are  not  in  accordance,  we  could  imagine  to  localize  the  spatial  region(s)  where 
the  reconstruction  is  no  longer  accurate.  The  domain  could  then  be  split  into  two  parts: 

•  a  subdomain  with  small  Kolmogorov  n-width  where  the  reconstruction  by  GEIM  is  still  accurate  enough. 

•  a  subdomain  with  big  Kolmogorov  n-width  where  the  accident  is  located  and  GEIM  is  no  longer  accurate. 
The  domain  is  computed  by  traditional  discretization  techniques  such  as  finite  elements  complemented  with 
Dirichlet  boundary  conditions  from  the  GEIM  reconstruction. 

Under  the  hypothesis  that  the  subdomain  affected  by  the  accident  is  small,  the  reconstruction  could  still  be  performed 
relatively  quickly,  preserving  the  real-time  aspect  of  our  device.  The  feasibility  of  decomposing  the  domain  and 
coupling  GEIM  with  other  approximations  has  been  explored  in  [1]  in  a  simple  Laplace  problem. 

Last  but  not  least,  it  would  also  be  interesting  to  explore  the  robustness  of  the  method  in  cases  where  one  or  several 
sensors  involved  in  the  GEIM  reconstmction  fail  either  by  being  mute  or  providing  a  wrong  answer. 
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Appendix  A. 

Proposition  Appendix  A.l.  Let  X  be  a  Hilbert  space  and  E,  F  two  subspaces  of  X.  Then,  Pe,f  —  where: 


Pe,f=  inf  sup,(e,/) 

eeE 

"^'1=' ll/ll 

/3f\e^  =  inf  sup(e,/). 

eeE^ 

11/11=1  IHI=1 


(A.l) 

(A.2) 


Proof.  Given  e  e  A  of  norm  unity,  we  introduce  f*  as 


/;  =  arg  sup(e,g). 

geF 

\\g\\=^ 

We  can  then  show  from  optimality  that  (e,  h)  -0  for  all  hm{q  e  F  \  (q,  f*)  -  0)  and  hence 

e  -  Af:  +  s  (A.3) 

for  some  A  eTl  and  e  e  F^  such  that  A^  +  ||e|p  =  1  (from  our  normalization  and  orthogonality).  We  then  deduce  from 
(A.3),  orthogonality,  and  Cauchy-Schwarz  that 

sup  (e,p)  -  A 

peF 

Ibihl 


and 


Hence, 


sup(e,p)  =  Hell. 

peF^ 


(  sup  (e,  p)f  +  { sup  (e,  p)f 


peF 

\\p\\=i 

thanks  to  our  normalization. 

We  may  now  note  from  (A.l)  and  (A.4)  that 

Pe,f  - 


peF^ 


inf  1  -  ( sup  (e,  p)  f 


1  -  ( sup  sup  (e,  p)  f 


eeE  peE 
M=l 


/I  -  ( sup  sup  (e,  p)  y- 

peF^  eeE 
Li=i  lkll=l 


=  1  I 

as  we  can  exchange  the  two  supremizer  operations. 

Finally,  we  dehne  a  second  inf-sup  constant, 

/3f^,e^  =  inf  sup(e,/). 

HF"-  eeE^ 

11/11=1  IHI=1 

We  can  repeat  the  procedure  above  —  E  goes  to  F-'-  and  F  goes  to  E-'-  —  to  hnd 

Pe^,e^  =  /l  -  (  sup  sup  (e,  p)  )2, 

-i  I  peF^  ee(Ey^ 

V  llpll=l  IHI=1 


and  hence  conclude  from  (A. 5)  and  (A. 7)  that 
since  (F''')'''  =  E. 


Pe,f  -  Pe^,e^ 


(A.4) 


(A.5) 


(A.6) 


(A.7) 
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Appendix  B. 


We  propose  here  a  practical  method  for  the  computation  of 


\\jMm\x 

ipeX  \m\x 


(B.l) 


The  strategy  consists  in  using  a  finite  element  Galerkin  projection  as  an  approximation  of  the  elements  of  X.  We 
therefore  propose  to  compute 

WSf M['/’]|lvf 

— nii - 

peVl  m\vt 

"  b 

as  a  surrogate  of  (B.l),  where  is  the  classical  continuous  finite  element  approximating  space  of  mesh  size  h  that 
involves  piece-wise  polynomials.  Let  S  =  spanj/ji, . . . ,  bf^]  be  a  basis  of  Vj  and  let  M  be  the  N  x  N  mass  matrix 
of  entries  Mjj  -  (bi,  bj)x,  1  <  i,  j  <  N.  For  any  ip  e  let 


(B.2) 

be  the  vector  of  coordinates  of  ip  in  the  basis  S.  In  coherence  with  these  notations,  for  any  I  <i  <  M,  the  vectors 

qi  =  (qij, qujf  and  Wi  =  (wi,;, . . . ,  (B.3) 

will  respectively  denote  the  Galerkin  projections  onto  of  the  interpolating  basis  functions  qi  e  X  and  of  the  Riesz 
representation  of  the  /-th  linear  functional,  cr,.  Furthermore,  let  be  the  N  x  M  matrix  such  that 


and  let  be  the  M  x  Af  matrix  such  that: 

c"  =  (Tiibj)  =  {Wi, bj)x,  Vl<i<M,  l<j<N. 

Finally,  we  recall  that  is  the  M  x  M  matrix  defined  in  Section  1  whose  entries  are 

B^j  =  (Tiiqj)  =  (wi,  qj)x,  VI  <  i  <  M,  \  <  j  <M. 

An  approximation  of  the  entries  of  and  can  easily  be  computed  by  using  the  finite  element  Galerkin  projections 
of  the  involved  functions: 

VI  </<M,  1  <;■<  Af 
Wi^Mqj,  VI  <i  <  M,  I  <  j  <  M. 

With  these  notations,  we  can  easily  prove 

Lemma  Appendix  B.l.  Let  T  be  the  N  x  N  symmetric  positive  definite  matrix: 

and  let  /Imax(T)  be  the  largest  eigenvalue  of  the  generalized  eigenvalue  problem 

[pind  (A,  x)  e  M  X  such  that: 

j  B.4) 

Tx  =  AMx. 


Then: 


IIJ'mMIU  rr 

max -  =  vdn 

pevl  II^IU  ^ 


(B.5) 
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Proof.  For  any  and  any  \  <  i  <  M: 


N 

o-ii^p)  =  f^(pjcri{bj)  =  ei^c’^(p, 
j=i 


where  et  is  the  i-ih  canonical  vector  of  dimension  M.  Furthermore,  if 

M 

JmIv]  =  ^  af{ip)qi 

i=\ 

is  the  generalized  interpolant  of  ip  in  dimension  M,  we  have: 


cTi  (Jum  =  ei^B^a,  VI  <  /  <  M, 

where  or  =  . . . ,  .  From  the  interpolation  property  stated  in  (3),  it  follows  that 

a  =  C^(p. 

Then,  the  finite  element  Galerkin  projection  of  the  interpolant  of  (B.6)  can  be  expressed  as: 

j'mM  ^  =  e"  cv 

Hence, 


max 

<fevl 


IIJ'mMIU 

II^IU 


'qM(bM)-Icm' 

T 

M 

<p 

\ 

<p^M<p 

> 

V'^max(T). 


(B.6) 


□ 

Remark  Appendix  B.2.  The  computation  ofA^ax  can  easily  be  performed  by,  e.g.,  the  power  method  scheme  applied 
to  the  matrix  T.  However,  note  that  the  evaluation  of  Am  with  formula  (B.5)  requires  the  construction  ofT,  which 
is  a  large  dense  matrix  of  dimension  N  X  N.  In  cases  where  the  storage  of  T  is  no  longer  possible,  the  Lebesgue 
constant  can  still  be  computed  with  formula  ( 30),  whose  evaluation  requires  the  construction  of  a  much  smaller  matrix 
of  dimension  M  X  M. 
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